home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / MacII / Accuracy / Accuracy.p next >
Encoding:
Text File  |  1988-05-23  |  16.3 KB  |  542 lines  |  [TEXT/ttxt]
  1. { ACCURACY.PAS: 8087 math accuracy test, v 1.7, Oct 1987
  2.     (c) 1987, PC Tech Journal and Ziff Communications Co.
  3.                     written by Jim Roberts.
  4.  
  5. Uncomment the correct constant declarations for your compiler.
  6. Insert correct version number in `compil' variable,
  7. name of hardware system & math coprocessor in `machin' variable.
  8.  
  9. }
  10. {$MC68020+}
  11. {$MC68881+}
  12. {$R+}
  13. {$U-}
  14.  
  15. PROGRAM Accuracy2;
  16.  
  17.   USES Memtypes, QuickDraw, OSIntf, ToolIntf, PackIntf, SANE;
  18.  
  19.   TYPE
  20.     accarray = array [1..5] of extended;
  21.  
  22.   CONST
  23.     minerr = 1.0E-17;
  24.     logmin = 17.0;
  25.     n = 10;
  26.     y = 1.0;
  27.     step = 0.2;
  28.     iter = 20;
  29.     itertrig = 5;
  30.     log10e = 0.43429448190325182765;
  31.     pi = 3.14159265358979323846;
  32.     piO2 = 1.57079632679489661923;
  33.     root2 = 1.4142135623730950488;
  34.     root3 = 1.7320508075688772935;
  35.     sqrtO2 = 0.7071067811865475244;
  36.  
  37.   VAR
  38.     i, j, k, l, m, ntest, z: integer;
  39.     a, b, c: array [1..n, 1..n] of extended;
  40.     sum, X: extended;
  41.     xx, zz, quot: extended;
  42.     a0, a1, d0, d1, frac: extended;
  43.     p, p2: extended;
  44.     th, err, logerr, diverr, val, funct: accarray;
  45.     testerr: array [1..10] of extended;
  46.     toterr: extended;
  47.     done: Boolean;
  48.     s: DecStr;
  49.     textWide: integer; { Maximum width of a character }
  50.     textHeight: integer; { Height of a character }
  51.     f: DecForm;
  52.  
  53.   Procedure Pause;
  54.   { Wait for mouse button click.  Also clear keyboard events. }
  55.     begin
  56.       while Button do SystemTask;
  57.       while not Button do SystemTask;
  58.       FlushEvents(keyDownMask + autoKeyMask + mDownMask, 0)
  59.     end; { Pause }
  60.  
  61.   Procedure DrawAt(X, y: integer;
  62.                    s: Str255);
  63.   { Draw string at this coordinate, similar to an x,y location }
  64.   { on a conventional computer terminal. }
  65.     const
  66.       Xoffset = 5; { Blank pixels in left border }
  67.       Yoffset = 12; { Blank pixels in top border }
  68.     begin
  69.       MoveTo(Xoffset + X * textWide, Yoffset + y * textHeight);
  70.       DrawString(s)
  71.     end; { DrawAt }
  72.  
  73.   Procedure NewFont(fontNumber, pointSize: integer);
  74.   { Select new font and adjust variables for DrawAt }
  75.     begin
  76.       TextFont(fontNumber); { Use new font }
  77.       TextSize(pointSize); { Use new size }
  78.       textHeight := pointSize + 2; { Calc new character height }
  79.       textWide := CharWidth('M') { Calc width of widest char }
  80.     end; { NewFont }
  81.  
  82.   Procedure filla;
  83.     var
  84.       i, j: integer;
  85.     begin
  86.       for i := 1 to n do
  87.         for j := 1 to n do
  88.           if i <> j then
  89.             a[i, j] := y
  90.           else
  91.             a[i, j] := X + y;
  92.     end;
  93.  
  94.   Procedure fillb;
  95.     var
  96.       i, j: integer;
  97.       f, d: extended;
  98.     begin
  99.       f := X + n * y;
  100.       d := 1.0 / (X * f);
  101.       for i := 1 to n do
  102.         for j := 1 to n do
  103.           if i <> j then
  104.             b[i, j] := - y * d
  105.           else
  106.             b[i, j] := ( - y + f) * d;
  107.     end;
  108.  
  109.   Procedure fillc;
  110.     var
  111.       i, j: integer;
  112.     begin
  113.       for i := 1 to n do for j := 1 to n do c[i, j] := 0.0;
  114.     end;
  115.  
  116.   Procedure mult;
  117.     var
  118.       i, j, k: integer;
  119.     begin
  120.       for i := 1 to n do
  121.         for j := 1 to n do
  122.           begin
  123.             sum := 0.0;
  124.             for k := 1 to n do sum := sum + a[i, k] * b[k, j];
  125.             c[i, j] := sum;
  126.           end;
  127.     end;
  128.  
  129.   Procedure sumit;
  130.     var
  131.       i, j: integer;
  132.     begin
  133.       sum := 0.0;
  134.       for i := 1 to n do c[i, i] := c[i, i] - 1.0;
  135.       for i := 1 to n do for j := 1 to n do sum := sum + abs(c[i, j]);
  136.     end;
  137.  
  138.   function osgn(n: integer): extended; {negative if n odd}
  139.     begin
  140.       if n = 2 * (n div 2) then
  141.         osgn := 1.0
  142.       else
  143.         osgn := - 1.0;
  144.     end;
  145.  
  146.   Procedure arith;
  147.   {TEST 1: well-conditioned combinatorial matrix times its inverse.}
  148.     var
  149.       i, k, l, m: integer;
  150.     begin
  151.       zz := 0.30; {factor used to control decrease of condition of matrix}
  152.       for l := 1 to 5 do
  153.         begin
  154.           xx := zz * (3 - l);
  155.           X := exp(xx / log10e); {slowly decreases conditioning}
  156.           filla;
  157.           fillb;
  158.           fillc;
  159.           mult;
  160.           sumit;
  161.           err[l] := sum / sqr(n); {error is average absolute error per element}
  162.           if err[l] > minerr then
  163.             logerr[l] := - ln(err[l]) * log10e
  164.           else
  165.             logerr[l] := logmin;
  166.           testerr[1] := testerr[1] + (logmin - logerr[l]);
  167.         end;
  168.       testerr[1] := testerr[1] / 5.0;
  169.       DrawAt(5, 21, '#1a: 10x10 matrix       ');
  170.       f.Style := FixedDecimal;
  171.       f.Digits := 1;
  172.       z := 25;
  173.       for i := 1 to 5 do
  174.         begin
  175.           z := z + 5;
  176.           Num2Str(f, logerr[i], s);
  177.           DrawAt(z, 21, s);
  178.         end;
  179.       f.Digits := 2;
  180.       Num2Str(f, testerr[1], s);
  181.       DrawAt(59, 21, s);
  182.       { infinite product for 1-x: run in reverse to test division }
  183.       sum := 0.0;
  184.       for l := 1 to 5 do
  185.         begin
  186.           xx := (1 - l) / 4.0;
  187.           zz := exp((xx - 2.0) / log10e); {increases number of factors for convergence}
  188.           xx := 1.0 - zz; { zz ≈ .01 => loss of 2 sig figs }
  189.    {The following formula for the number of factors is designed to give
  190.       sufficient accuracy, while avoiding underflow in the powers of xx.
  191.       It gives a more uniform computation from compiler to compiler.}
  192.           m := 12 + l;
  193.           quot := 1.0;
  194.           for k := 1 to m do
  195.             begin
  196.               quot := quot / (1.0 + xx);
  197.               xx := xx * xx;
  198.             end;
  199.           err[l] := abs(1.0 - quot / zz) * 0.01;
  200.           { factor of 0.01 to compensate for cancellation errors }
  201.           if err[l] > minerr then
  202.             diverr[l] := - ln(err[l]) * log10e
  203.           else
  204.             diverr[l] := logmin;
  205.           sum := sum + (logmin - diverr[l]);
  206.           logerr[l] := diverr[l]; {needed for later average}
  207.         end;
  208.       sum := sum / 5.0;
  209.       DrawAt(5, 22, '#1 : infinite product   ');
  210.       f.Digits := 1;
  211.       z := 25;
  212.       for i := 1 to 5 do
  213.         begin
  214.           z := z + 5;
  215.           Num2Str(f, diverr[i], s);
  216.           DrawAt(z, 22, s);
  217.         end;
  218.       f.Digits := 2;
  219.       Num2Str(f, sum, s);
  220.       DrawAt(59, 22, s);
  221. { test continued fraction for tangent against exact values for five angles:
  222.    this is a test of division and subtraction, not of the tangent.}
  223.       th[1] := pi / 12.0;
  224.       th[2] := pi / 6.0;
  225.       th[3] := pi / 4.0;
  226.       th[4] := pi / 3.0;
  227.       th[5] := 5.0 * pi / 12.0;
  228.       val[1] := 2.0 - root3;
  229.       val[2] := 1.0 / root3;
  230.       val[3] := 1.00;
  231.       val[4] := root3;
  232.       val[5] := 2.0 + root3;
  233.       sum := 0.0;
  234.       m := 8; { this number of iterations gives sufficient accuracy }
  235.       for l := 1 to 5 do
  236.         begin
  237.           a0 := 2.0 * m + 1.0;
  238.           p2 := th[l];
  239.           p := sqr(p2);
  240.           d0 := a0 - p / (a0 + 2.0);
  241.           for k := 1 to m do
  242.             begin
  243.               a1 := a0 - 2.0;
  244.               d1 := a1 - p / d0;
  245.               a0 := a1;
  246.               d0 := d1;
  247.             end;
  248.           frac := p2 / d0;
  249.           funct[l] := frac;
  250.         end;
  251.       for l := 1 to 5 do
  252.         begin
  253.           err[l] := abs(1.0 - val[l] / funct[l]);
  254.           if err[l] > minerr then
  255.             diverr[l] := - ln(err[l]) * log10e
  256.           else
  257.             diverr[l] := logmin;
  258.           sum := sum + (logmin - diverr[l]);
  259.         end;
  260.       sum := sum / 5.0;
  261.       DrawAt(5, 23, '#1 : continued fraction ');
  262.       f.Digits := 1;
  263.       z := 25;
  264.       for i := 1 to 5 do
  265.         begin
  266.           z := z + 5;
  267.           Num2Str(f, diverr[i], s);
  268.           DrawAt(z, 23, s);
  269.         end;
  270.       f.Digits := 2;
  271.       Num2Str(f, sum, s);
  272.       DrawAt(59, 23, s);
  273.       for i := 1 to 5 do
  274.         begin
  275.           logerr[i] := 0.5 * (logerr[i] + diverr[i]);
  276.           testerr[2] := testerr[2] + (logmin - logerr[i]);
  277.         end;
  278.       testerr[2] := testerr[2] / 5.0;
  279.       DrawAt(5, 24, '#1b: division average   ');
  280.       f.Digits := 1;
  281.       z := 25;
  282.       for i := 1 to 5 do
  283.         begin
  284.           z := z + 5;
  285.           Num2Str(f, logerr[i], s);
  286.           DrawAt(z, 24, s);
  287.         end;
  288.       f.Digits := 2;
  289.       Num2Str(f, testerr[2], s);
  290.       DrawAt(59, 24, s);
  291.     end;
  292.  
  293.   Procedure trig; {TEST 2: first, errors in some sine identities }
  294.     var
  295.       i, j, k, l: integer;
  296.     begin
  297.       for l := 1 to 5 do logerr[l] := 0.0;
  298.       for j := 1 to itertrig do
  299.         begin
  300.           p := j - 1;
  301.           th[1] := pi / 12.0 + p * pi;
  302.           th[2] := pi / 6.0 + p * pi;
  303.           th[3] := pi / 4.0 + p * pi;
  304.           th[4] := pi / 3.0 + p * pi;
  305.           th[5] := 5.0 * pi / 12.0 + p * pi;
  306.           val[1] := osgn(j - 1) * root2 * (root3 - 1.0) * 0.25;
  307.           val[2] := osgn(j - 1) * 0.5;
  308.           val[3] := osgn(j - 1) * sqrtO2;
  309.           val[4] := osgn(j - 1) * 0.5 * root3;
  310.           val[5] := osgn(j - 1) * root2 * (root3 + 1.0) * 0.25;
  311.           for l := 1 to 5 do funct[l] := sin(th[l]);
  312.           for l := 1 to 5 do
  313.             begin
  314.               err[l] := abs(1.0 - val[l] / funct[l]);
  315.               if err[l] > minerr then
  316.                 logerr[l] := logerr[l] - ln(err[l]) * log10e
  317.               else
  318.                 logerr[l] := logerr[l] + logmin;
  319.             end;
  320.         end;
  321.       for l := 1 to 5 do logerr[l] := logerr[l] / itertrig;
  322.       for l := 1 to 5 do testerr[3] := testerr[3] + (logmin - logerr[l]);
  323.       testerr[3] := testerr[3] / 5.0;
  324.       DrawAt(5, 25, '#2a: sine function      ');
  325.       f.Digits := 1;
  326.       z := 25;
  327.       for i := 1 to 5 do
  328.         begin
  329.           z := z + 5;
  330.           Num2Str(f, logerr[i], s);
  331.           DrawAt(z, 25, s);
  332.         end;
  333.       f.Digits := 2;
  334.       Num2Str(f, testerr[3], s);
  335.       DrawAt(59, 25, s);
  336.  
  337.       { compare sin()/cos() with exact values for tan() }
  338.       for l := 1 to 5 do logerr[l] := 0.0;
  339.       for j := 1 to itertrig do
  340.         begin
  341.           p := j - 1;
  342.           th[1] := pi / 12.0 + (j - 1) * pi;
  343.           th[2] := pi / 6.0 + (j - 1) * pi;
  344.           th[3] := pi / 4.0 + (j - 1) * pi;
  345.           th[4] := pi / 3.0 + (j - 1) * pi;
  346.           th[5] := 5.0 * pi / 12.0 + (j - 1) * pi;
  347.           val[1] := 2.0 - root3;
  348.           val[2] := 1.0 / root3;
  349.           val[3] := 1.00;
  350.           val[4] := root3;
  351.           val[5] := 2.0 + root3;
  352.           for l := 1 to 5 do funct[l] := sin(th[l]) / cos(th[l]);
  353.           for l := 1 to 5 do
  354.             begin
  355.               err[l] := abs(1.0 - val[l] / funct[l]);
  356.               if err[l] > minerr then
  357.                 logerr[l] := logerr[l] - ln(err[l]) * log10e
  358.               else
  359.                 logerr[l] := logerr[l] + logmin;
  360.             end;
  361.         end;
  362.       for l := 1 to 5 do logerr[l] := logerr[l] / itertrig;
  363.       for l := 1 to 5 do testerr[4] := testerr[4] + (logmin - logerr[l]);
  364.       testerr[4] := testerr[4] / 5.0;
  365.       DrawAt(5, 26, '#2b: tangent|sin()/cos()');
  366.       f.Digits := 1;
  367.       z := 25;
  368.       for i := 1 to 5 do
  369.         begin
  370.           z := z + 5;
  371.           Num2Str(f, logerr[i], s);
  372.           DrawAt(z, 26, s);
  373.         end;
  374.       f.Digits := 2;
  375.       Num2Str(f, testerr[4], s);
  376.       DrawAt(59, 26, s);
  377.             
  378.       { compare arctan() with tan() for consistency }
  379.       for l := 1 to 5 do logerr[l] := 0.0;
  380.       for j := 1 to iter do
  381.         begin
  382.           for l := 1 to 5 do th[l] := (5 * j + l - 5) * piO2 / (5 * iter + 1);
  383.           for l := 1 to 5 do val[l] := sin(th[l]) / cos(th[l]);
  384.           for l := 1 to 5 do funct[l] := arctan(val[l]);
  385.           for l := 1 to 5 do
  386.             begin
  387.               err[l] := abs(1.0 - th[l] / funct[l]);
  388.               if err[l] > minerr then
  389.                 logerr[l] := logerr[l] - ln(err[l]) * log10e
  390.               else
  391.                 logerr[l] := logerr[l] + logmin;
  392.             end;
  393.         end;
  394.       for l := 1 to 5 do logerr[l] := logerr[l] / iter;
  395.       for l := 1 to 5 do testerr[5] := testerr[5] + (logmin - logerr[l]);
  396.       testerr[5] := testerr[5] / 5.0;
  397.       DrawAt(5, 27, '#2c: arctan function    ');
  398.       f.Digits := 1;
  399.       z := 25;
  400.       for i := 1 to 5 do
  401.         begin
  402.           z := z + 5;
  403.           Num2Str(f, logerr[i], s);
  404.           DrawAt(z, 27, s);
  405.         end;
  406.       f.Digits := 2;
  407.       Num2Str(f, testerr[5], s);
  408.       DrawAt(59, 27, s);
  409.     end;
  410.  
  411.   Procedure transc; {TEST 3: ln() and exp() for consistency}
  412.     var
  413.       i, j, l: integer;
  414.     begin
  415.       for l := 1 to 5 do logerr[l] := 0.0;
  416.       for j := 1 to iter do
  417.         begin
  418.           for l := 1 to 5 do th[l] := (5 * j + l - 5) * step;
  419.           for l := 1 to 5 do val[l] := exp(th[l]);
  420.           for l := 1 to 5 do funct[l] := ln(val[l]);
  421.           for l := 1 to 5 do
  422.             begin
  423.               err[l] := abs(1.0 - th[l] / funct[l]);
  424.               if err[l] > minerr then
  425.                 logerr[l] := logerr[l] - ln(err[l]) * log10e
  426.               else
  427.                 logerr[l] := logerr[l] + logmin;
  428.             end;
  429.         end;
  430.       for l := 1 to 5 do logerr[l] := logerr[l] / iter;
  431.       for l := 1 to 5 do testerr[6] := testerr[6] + (logmin - logerr[l]);
  432.       testerr[6] := testerr[6] / 5.0;
  433.       DrawAt(5, 28, '#3a: ln() & exp()       ');
  434.       f.Digits := 1;
  435.       z := 25;
  436.       for i := 1 to 5 do
  437.         begin
  438.           z := z + 5;
  439.           Num2Str(f, logerr[i], s);
  440.           DrawAt(z, 28, s);
  441.         end;
  442.       f.Digits := 2;
  443.       Num2Str(f, testerr[6], s);
  444.       DrawAt(59, 28, s);
  445.     end;
  446.  
  447.   Procedure roots; { sqrt() and sqr() identities }
  448.     var
  449.       i, j, l: integer;
  450.     begin
  451.       for l := 1 to 5 do logerr[l] := 0.0;
  452.       for j := 1 to iter do
  453.         begin
  454.           for l := 1 to 5 do th[l] := (5 * j + l - 5) * step;
  455.           for l := 1 to 5 do val[l] := sqrt(th[l]);
  456.           for l := 1 to 5 do funct[l] := sqr(val[l]);
  457.           for l := 1 to 5 do
  458.             begin
  459.               err[l] := abs(1.0 - th[l] / funct[l]);
  460.               if err[l] > minerr then
  461.                 logerr[l] := logerr[l] - ln(err[l]) * log10e
  462.               else
  463.                 logerr[l] := logerr[l] + logmin;
  464.             end;
  465.         end;
  466.       for l := 1 to 5 do logerr[l] := logerr[l] / iter;
  467.       for l := 1 to 5 do testerr[7] := testerr[7] + (logmin - logerr[l]);
  468.       testerr[7] := testerr[7] / 5.0;
  469.       DrawAt(5, 29, '#3b: sqrt() & sqr()      ');
  470.       f.Digits := 1;
  471.       z := 25;
  472.       for i := 1 to 5 do
  473.         begin
  474.           z := z + 5;
  475.           Num2Str(f, logerr[i], s);
  476.           DrawAt(z, 29, s);
  477.         end;
  478.       f.Digits := 2;
  479.       Num2Str(f, testerr[7], s);
  480.       DrawAt(59, 29, s);
  481.     end;
  482.  
  483.   Procedure Initialize;
  484.   { Perform one-time setup work. }
  485.     var
  486.       r: Rect;
  487.       MyNewWindow: WindowPtr;
  488.     begin
  489.       InitGraf(@thePort); { initialize QuickDraw }
  490.       InitFonts; { " Font Manager }
  491.       InitWindows; { " Window Manager }
  492.       InitCursor; { make cursor an arrow }
  493.       done := False;
  494.       SetRect(r, 30, 60, 600, 400);
  495.       MyNewWindow := NewWindow(nil, r, 'Accuracy Tester', True, documentProc, Pointer( - 1), True,
  496.                                0);
  497.       MyNewWindow^.txFace := [bold];
  498.       MyNewWindow^.txFont := SystemFont;
  499.       windowPeek(MyNewWindow)^.refCon := Longint(NewString('Contents of window 1'));
  500.     end;
  501.  
  502.   Procedure Header;
  503.     begin
  504.       NewFont(monaco, 9);
  505.       TextFace([]);
  506.       DrawAt(5, 6, 'ACCURACY: extended tester: MPW Pascal: Macintosh II');
  507.       DrawAt(5, 7, '  V 1.7 (c) 1987, PC Tech Journal and Ziff Communications Co.');
  508.       DrawAt(5, 8, '                  written by Jim Roberts.');
  509.             DrawAt(5,9,'    Modified for Macintosh & Turbo Pascal by Morry Hodges, CIS 74766,3426');
  510.       DrawAt(5, 10, 'Test 1 checks multiplication and addition, then division and subtraction.');
  511.       DrawAt(5, 11,
  512.              'Test 2 measures the accuracy of the Turbo trig functions (there is no tan()).');
  513.       DrawAt(5, 12, 'Test 3 finds the truncation error in some exponential and sqrt identities.');
  514.       DrawAt(5, 13, 'ACCURACY is the rounded negative log of error.  Program may exit abnormally.');
  515.       DrawAt(5, 14, 'NOTE: an increase of 1 in the rating means a factor of TEN less accurate.');
  516.       DrawAt(5, 16, 'Interpretation  <0.0 - 0.5 => Excellent     1.0 - 1.5 => Fair');
  517.       DrawAt(5, 17, '  of RATING:     0.5 - 1.0 => Good          1.5 <     => Poor');
  518.       DrawAt(5, 18, '');
  519.       DrawAt(5, 19, '      TESTS                       ACCURACY           RATING ');
  520.     end;
  521.  
  522.   begin { program ACCURACY begins here }
  523.     repeat
  524.       Initialize;
  525.       Header;
  526.       for i := 1 to 10 do testerr[i] := 0.0;
  527.       arith;
  528.       trig;
  529.       transc;
  530.       roots;
  531.       ntest := 7;
  532.       toterr := 0.0;
  533.       for i := 1 to ntest do toterr := toterr + testerr[i];
  534.       toterr := toterr / ntest;
  535.       DrawAt(5, 31, 'Overall rating:');
  536.       Num2Str(f, toterr, s);
  537.       DrawAt(23, 31, s);
  538.       DrawAt(5, 33, 'Click mouse to finish');
  539.       Pause;
  540.     until Button;
  541.   end.
  542.